#include "Cf252.hh"
#include "G4ParticleTable.hh"
#include "G4Neutron.hh"
#include "G4Gamma.hh"


Cf252::Cf252(G4ThreeVector p)
: neutronSource(0)
{
	position = p;
	createFissionSource(); 
	createAlphaSource(); 
	
}

Cf252::~Cf252()
{
	//nothing extra to do.
}

void Cf252::GeneratePrimaryVertex(G4Event* evt)
{
	G4int numberOfFissionNeutrons = sampleNeutronNumberDist();
	neutronSource->SetNumberOfParticles(numberOfFissionNeutrons);
	BranchingPrimaryGenerator::GeneratePrimaryVertex(evt);
}

void Cf252::createAlphaSource()
{
	
	/* insert code here. 
	One way to do this would be to create a BranchingPrimaryGenerator 
	"alphaBPG" and fill it with instances of G4GeneralParticleSource. 
	You can set G4GeneralParticle source to isotropically emit an alpha 
	its corresponding gammas. The branching ratios for each alpha+gammas
	set can be set when you add them to "alphaBPG".

	
	Finally add the BranchingPrimaryGenerator "alphaBPG" you just created to 
	model the alpha decay of Cf252 to the generatorList here with the method
	this->addPrimaryGenerator(alphaBPG,branchingRatio);

	Use the position vector in this class to set the position of the sources.
	
	*/
	
}

void Cf252::createFissionSource()
{
//Creating a general particle source
G4GeneralParticleSource* fissionSource = new G4GeneralParticleSource();

// set the fissionSource to generate multiple vertexes, so it will generate the
// neutrons and the prompt fission gammas instead of just one or the other.
fissionSource->SetMultipleVertex(true);

//add fission neutron info to fission source;
configureFissionSourceForNeutrons(fissionSource);

//add prompt gamma info...~10.3 gammas/fission
configureFissionSourceForGammas(fissionSource);
//add fission source to generator list for Cf252 (from BranchingPrimaryGenerator)
this->addPrimaryGenerator(fissionSource, //source 
						  1.0);//1-0.969); //branching ratio for spontaneous fission.


//

	/* insert code here
	One way to do this would be to create a G4GeneralParticleSource and
	create two of G4SingleParticleSource in it. You can configure these sources
	to emit isotropically and you can turn on the SetMultipleVertex option
	to cause both the sources you add to be generated in an event. These 
	two sources would correspond to the prompt gammas and the fission neutrons.
	G4SingleParticleSource lets you set an energy distribution.

	Finally add the G4GeneralParticleSource to model the spontaneous fision of
	Cf252 to the generatorList with the method 
	this->addPrimaryGenerator(fissionGPS,branchingRatio);

	Use the position vector in this class to set the position of the sources.

	*/
}



void Cf252::configureFissionSourceForNeutrons(G4GeneralParticleSource* fissionSource)
{
	//creating a neutron source
	fissionSource->AddaSource(1.0);
	// get pointer to neutron single particle source
	this->neutronSource = fissionSource->GetCurrentSource();


	//set the particle definition
	//G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable() ->FindParticle("neutron");
	G4ParticleDefinition* neutronDef = G4Neutron::Definition();
	neutronSource->SetParticleDefinition(neutronDef);
	
	//This sets the source to isotropic. default is planar.
	neutronSource->GetAngDist()->SetAngDistType("iso");


	/*** the code below is for setting   ////
	//// the neutron energy distribution ***/

	G4SPSEneDistribution* eDistGen = 0;
	eDistGen = neutronSource->GetEneDist(); //returns pointer to G4SPSEneDistribution
	
	//tell energy distribution we want to use our own energy spectrum
	eDistGen->SetEnergyDisType("User");
	//set the random generator. Don't know why it doesn't create one automatically.
	//eDistGen->SetBiasRndm(new G4SPSRandomGenerator());
	//choose neutron energy distribution from Watt dist from Radiation Shielding textbook
	double a, b, c, chi, E, Emax, dE;
	a = 0.64;
	b = 1.175;
	c = 1.0401;
	Emax = 20;
	dE = 0.1;
	E = 0.05;

	while(E < Emax)
	{
		chi = a*exp(-E/b)*sinh(sqrt(c*E));
		//G4cout << "chi = " << chi << "\n";
		//giving the histogram values (don't know if this is normallized yet)
		eDistGen->UserEnergyHisto(G4ThreeVector(E*MeV,chi,0));
		E = E + dE;
	}
    
	//eDistGen->UserEnergyHisto(G4ThreeVector(10*MeV,1.,0));
	//eDistGen->UserEnergyHisto(G4ThreeVector(10*MeV,0.5,0));

	//*** done setting the neutron energy distribution ***

	//set the initial number of fission neutrons.
	G4int initialNumFissionNeuts = sampleNeutronNumberDist();
	neutronSource->SetNumberOfParticles( initialNumFissionNeuts );
	 
	
	
}



G4int Cf252::sampleNeutronNumberDist()
{
	// Number of fission neutrons is adequately 
	// described by a gaussian. These parameters are from
	// "Distributions of Fission Neutron Numbers" by
	// James Terrell.
	G4double mean, std, gaussianSample, decimalPart;
	G4int integerPart;
	
	//sample from a gaussian distribution with the parameters
	//from Terrell.
	gaussianSample = G4RandGauss::shoot(mean = 3.86, std = 1.207);
	
	//just some sanity numbers since we are using a Gaussian
	if(gaussianSample <=0)
	{
		return 0;
	}
	else if(gaussianSample > 8)
	{
		return 8;
	}
	//get the integer part of the floating point number;
	integerPart = (int)gaussianSample;

	//check to see if we need to round up.
	decimalPart = gaussianSample - integerPart;
	if(decimalPart >= 0.5)
	{
		integerPart++;
	}
	
	//return this number as the number of fission neutrons.
	return integerPart ;
}

void Cf252::configureFissionSourceForGammas(G4GeneralParticleSource* fissionSource)
{
	//Adds a New Single Particle Source to Fission Source
	fissionSource->AddaSource(1.0);
	//Getting a pointer from fission source to the new single particle source
	G4SingleParticleSource* gammaSource = fissionSource->GetCurrentSource();
	//Telling the newly created single particle source ""gammaSource" to emmit gamma rays
	G4ParticleDefinition* gammaDef = G4Gamma::Definition();
	gammaSource->SetParticleDefinition(gammaDef);
	
	//This sets the source to isotropic. default is planar.
	gammaSource->GetAngDist()->SetAngDistType("iso");

	G4SPSEneDistribution* eDistGen = gammaSource->GetEneDist();
	
	//tell eDistGen we want to use our own energy distribution
	eDistGen->SetEnergyDisType("User");
	//set the random generator. don't know why it's not automatic.
	//eDistGen->SetBiasRndm(new G4SPSRandomGenerator());
	//choose gamma energy distribution from U235-promp fission
	//from Radiation Shielding textbook pg 95
        double N, E, Emax, dE;
        Emax = 10.5;
	//dE should be .0001 or smaller in order to get as near as possible to a normalized pdf
        dE = 0.1;
        E = 0.1;

        while(E < Emax)
        {
		if (E <= 0.6)
		{
			N=6.6;///(7.7022/dE);
		}
		else if (E <= 1.5)
        	{
			N=20.2*exp(-1.78*E);///(7.7022/dE);
		}
    	else if (E <= 10.5)
        {
			N=7.2*exp(-1.09*E);///(7.7022/dE);
		}
		//G4cout << E << " = E   " << N << " = N\n";
		eDistGen->UserEnergyHisto(G4ThreeVector(E*MeV,N,0));
		E = E + dE;
        }

	//eDistGen->UserEnergyHisto(G4ThreeVector(5*MeV,1.,0));
	
	
	
	
	
	
	gammaSource->SetNumberOfParticles(10);
	
	

}
